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INTRODUCTION 


Obi ective.  The  specific  objective  is  to  formulate  a  viscoplastic  consti¬ 
tutive  model  which  incorporates  the  so-called  "CAP75"  plasticity  model  for 
geological  materials  (1).  Included  in  this  objective  is  the  development  of  a 
numerical  solution  algorithm  (computer  program)  which  "exercises"  the  visco¬ 
plastic  model  for  general  states  of  strain  loading  (or  stress  loading)  time 
histories. 

Background .  Geological  constitutive  models  based  on  inviscid  plasticity 
concepts,  such  as  CAP75,  are  inherently  independent  of  real  time.  Experimental 
evidence,  however,  indicates  a  strong  time-dependent  behavior  for  most  soils 
and  rocks  (2, 3,4, 5),  thus  motivating  the  above  objective.  The  coupling  of 
plastic  behavior  with  time  dependency  comes  under  the  general  heading  of  vis- 
coolasticity.  Like  religion,  viscoplasticity  has  a  variety  of  "sects"  each 
with  its  own  fervent  disciples,  dogmas,  claims  and  counterclaims.  Some  per¬ 
spective  on  the  variety  of  viscoplastic  formulations  can  be  gained  by  contrast¬ 
ing  the  endochronic  approach  of  Valinis  (6)  to  the  elastic/viscopiastic  approach 
of  Perzyna  (7).  The  genesis  of  the  former  may  be  envisioned  as  a  modification 
to  classical  viscoelasticity  wherein  plastic-like  behavior  is  introduced  by 
means  of  "intrinsic  time",  a  time  scale  considered  to  be  a  property  of  the  material. 
No  plastic  yield  function  is  introduced.  Alternatively,  Perzyna's  elastic/vis¬ 
copiastic  approach  is  a  modification  of  classical  plasticity  wherein  viscous-like 
behavior  is  introduced  by  a  time  rate  flow  rule  which  employs  a  plastic  yield 
function.  Other  viscoplastic  formulations  include  developments  by  Katona  (8), 

Bodner  and  Parton  (9) ,  and  Phillips  and  Wu  (10) . 

Any  viscoplastic  formulation  should  have  the  following  attributes: 

(1)  capability  of  simulating  observed  material  behavior  over  a  wide  range  of 
loading  conditions;  (2)  satisfaction  of  thermodynamic  restrictions;  (3)  feasible 
techniques  for  parameter  identification;  and  (4)  readily  adaptable  and  efficient 
in  numerical  solution  procedure  (e.g.  finite  element  methods).  With  regard  to 
the  last  attribute,  the  numerical  solution  strategy  is  clearly  of  utmost  impor¬ 
tance,  i.e.  it  makes  no  sense  to  utilize  sophisticated  constitutive  theories  if 
the  accuracy  of  the  solution  is  highly  questionable.  At  the  same  time,  however, 
some  compromise  in  accuracy  must  be  made  to  achieve  computational  efficiency. 


In  1974  Zienkiewlcz  and  Cormeau  (II)  presented  a  finite  element  algorithm 
for  a  viscoplastic  model  of  the  Perzyna  type  using  an  explicit  forward  difference 
scheme  for  time  integration.  Computationally  the  explicit  algorithm  is  simple, 
requiring  only  an  elastic  stiffness  matrix  to  be  assembled  and  triangularized 
at  the  outset.  Nonlinearities  are  accomodated  on  the  right-hand-side  at  the 
element  level  with  (if  desired)  a  mid-interval  iteration  and/or  next-interval 
equilibrium  correction  (12).  Unfortunately,  a  major  drawback  is  the  concern  for 
numerical  stability.  This  places  a  stringent  limit  on  the  size  of  the  time 
step  and,  thereby,  compromises  overall  efficiency.  Cormeau  (13)  presented  a  priori 
time  step  limits  to  avoid  instability  for  a  restricted  class  of  Perzyna-type 
viscoplastic  models.  However,  for  more  general  forms  (e.g.  viscoplastic  models 
with  hardening)  a  priori  time  step  limits  are  not  established. 

Numerical  stability  concerns  can  be  eliminated  by  employing  an  implicit 
time  integration  scheme.  Hughes  and  Taylor  (14)  showed  that  a  one  parameter 
(9)  Crank-Nicolson  integration  scheme  provides  unconditional  stability  for  0.5 
9  _<  1  when  applied  to  FEM  models  with  Perzyna-type  viscoplasticity.  For  a 
particular  problem,  they  demonstrated  accuracy  was  maintained  for  implicit 
time-step  sizes  10^  times  greater  than  the  critical  time-step  size  in  the 
explicit  algorithm. 

The  computational  penalty  for  implicit  integration  schemes  is  the  require¬ 
ment  to  solve  nonlinear  algebraic  equations  for  each  time  step.  Here  a  variety 
of  efficiency  vs.  accuracy  alternatives  are  possible  based  on  various  versions 
of  the  Newton-Raphson  technique  (15,16). 

Scope  and  Approach.  In  meeting  the  stated  objective,  a  Perzyna-type  visco¬ 
plastic  formulation  is  adopted.  Motivations  for  this  choice  are;  (1)  the  formu¬ 
lation  is  well  accepted  and  well  used,  and  (2)  the  incorporation  of  the  CAP75 
plasticity  model  (or  any  plasticity  model)  is  theoretically  straightforward. 

In  contrast,  the  endochronic  approach  to  viscoplasticity  continues  to  draw 
criticism  and  does  not  lend  itself  to  direct  incorporation  of  plasticity  models. 

For  this  study,  a  one  parameter  Crank-Nicolson  time  integration  scheme  is 
used  which  provides  options  for  explicit  or  implicit  algorithms.  For  the  im¬ 
plicit  algorithm,  nonlinear  algebraic  equations  are  solved  "exact-as-desired" 
by  an  iterative  Newton-Raphson  technique  with  special  modifications  for  plastic 
hardening . 
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For  clarity  and  insight,  the  presentation  begins  with  a  one-dimensional 
viscoplastic  formulation  along  with  some  exact  solutions  illustrating  the  in¬ 
fluence  of  viscoplastic  model  parameters.  The  numerical  time  integration  is 
also  introduced  at  the  one-dimensional  level  illustrating  the  nature  of  explicit 
and  implicit  algorithms.  Next,  the  general  multi-dimensional  viscoplastic  for¬ 
mulation  is  presented  together  with  the  numerical  solution  strategy.  Finally, 
the  formulation  is  specialized  for  the  CAP75  plasticity  model  and  example  re¬ 
sults  are  discussed. 

The  Appendix  contains  user  input  instructions  for  the  computer  program 
VP DRV R  which  is  a  viscoplastic  version  of  the  plasticity  program  CAPDRVR  (17) . 


ONE-DIMENSIONAL  VISCOPLASTICITY 

Although  one  dimensional  viscoplastic  stress-strain  relationships  have 
limited  practical  application,  they  provide  an  excellent  starting  point  for 
introducing  fundamental  theoretical  concepts,  behavior  characteristics,  and 
solution  strategies. 

Conceptual  Formulation.  Figure  1  illustrates  a  simple  viscoplastic 
rheological  model  composed  of  an  elastic  element  (e.g.  spring),  viscous  element 
(e.g.  dashpot),  and  a  plastic  element  (e.g.  slider).  In  such  a  rheological 
construction  we  first  identify  local  stress-strain  relations  of  the  component 
elements.  As  s  somewhat  general  example,  we  will  specify  the  following  forms. 

For  the  elastic  element  let: 


where  i  ,  c  *  stress  and  strain  in  elastic  element 
e  e 

E(c  )  *  modulus  function  of  elastic  strain, 

(linear  case  E(c  )  *  Ertc  ). 

e  0  e 

In  Equation  (1),  the  inverse  is  required  to  exist,  i.e.,  e  »  E~*(a  ). 

e  e 


I 
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The  viscous  stress-strain  relation  is  represented  in  power  form  as: 


/  -N 

/i°ji 


where 


Y  \“J  Sgn  (°v> 


0 

Y 


rate  of  viscous  straining,  dev/dt. 


viscous  stress. 

normalizing  constant  (units  of  stress), 
fluidity  parameter  (units  of  inverse  time). 


exponent  power  (N  >  0) . 


(2) 


sgn  (jv) 


.  +l,o  >0 
1  v 

-1,  a  <0 
v 


Note  the  linear  form  of  Equation  2  is  given  by  M  =  1  representing  a  simple 
dashpot  model.  In  such  a  case  the  dashpot  coefficient  is  u  *  jg/y,  (i.e. 

'-v  =  V  In  Che  m0re  §eneral  case  (N  *  1)  Equation  2  is  a  nonlinear  creep 

law  implying  chat  the  viscous  strain  rate  magnitude  is  proportional  to  some 
power  of  che  viscous  stress,  and  the  strain  rate  direction  (sign)  is  governed  by 
the  viscous  stress  sign. 


Lastly,  the  plastic  element  is  defined  by  means  of  an  isotropic  hardening 
plastic  yield  function  f; 


f  (a  ,  s  ) 

P  P 


ja  |  -  k(e  )  <  0 

P  P  “ 


(3) 


where 


7  ,  £ 

P  P 


stress  and  strain  in  plastic  element. 


k(e  ) 
P 


positive  valued  hardening  function  (for  non  hardening 

k(E  )  »  kn) . 
p  U 


In  conjunction  with  the  above,  we  adopt  the  classical  plasticity  assumptions  that 
plastic  flow  (increments  of  plastic  strain)  can  only  occur  when  the  equality  holds 
in  Equation  3,  i.e.; 


(4) 


0,  if  f(a  ,  e  )  <  0, 

P  P 

|de  |  sgn(a  ) ,  if  f(cr  ,  e  )  -  0, 

P  P  P  P 

and  -6gn(o  )do  >  0 
P  P 


Note,  f(o  ,  e  )  >  0  is  not  an  admissible  state  when  a  is  used  as  an  argument 
P  P  P 

in  f.  This  will  be  elaborated  subsequently. 

In  addition  to  the  above  component  relationships,  the  model  in  Figure  1 
inherently  implies  the  following  equilibrium  and  comparability  relationships 
which  must  always  hold  true. 


a  ■  a  »  a  +  a 
e  v  p 


(5) 


£ 


e 

e 


+  £ 

vp 


(6) 


C. 


vp 


£  * 

P 


c. 

V 


(7) 


sgn(a)  -  sgn  (a  )  -  sgn(a  ) 
P  v 


(8) 


where  a  *  total  stress 
£  *  total  strain 

*  viscoplastic  strain 

Ultimately  the  objective  is  to  utilize  Equations  1  through  8  to  get  a 
differential  equation  relating  total  stress  to  total  strain,  i.e.,  a  constitu¬ 
tive  relation.  To  this  end.  Equations  5  and  8  give,  |o  |  *  |<j|  -  jo^j,  which 
when  inserted  into  Equation  3  with  allows  the  left  hand  side  of  Equation 

3  to  be  replaced  by: 


f  (a 


£  ) 
vp 


f(a’  £vr,) 

vp 


(9) 

0  in  which  case 


By  virtue  of  Equation  U  plastic  flow  occurs  when  f(o  ,  t  ) 

P  VP 

Equation  9  gives  ja  |  "  f(a,  e  ).  Introducing  this  into  Equation  2  with 


6 


gives; 


-  c.  , 

v  vp 


e  «  y$(f )  sgn(a) 
vp 


where, 


when  f  >  0 


when  f  <  0 


and,  f(a,  £  ) 

vp 


k(e  ) 
vp 


Equation  10  is  a  fundamental  viscoplastic  flow  rule  of  the  Perzyna  type  (7)  which 
could  have  been  simply  stated  as  opposed  to  the  model  development  given  here. 
However,  certain  insights  from  the  model  development  are  useful.  For  example, 
the  so-called  "dynamic  yield  function",  given  by  Equation  12,  is  greater  than 
zero  during  plastic  flow  by  the  amount  ]a  *,  whereas  the  classical  plasticity 
yield  function  is  zero  during  plastic  flow  (Equation  3).  If  loading  is  held 
constant,  steady  state  conditions  are  eventually  reached,  =  £  =  0,  so  that, 

)  =  f(r,  s  )  -  0. 
p  vp  vp 

To  complete  the  constitutive  development  Equations  i,  5  and  6  are 
combined  to  give: 


Taking  the  time  derivative  of  Equation  13  and  replacing  £  by  Equation  10 
gives  me  final  constitutive  form: 


S(-;  -  vj> ( f )  sgn(j)) 


In  the  above,  f  contains  J  and  £  as  an  argument,  i.e.  f 

vp 


we  have  by  Equations  1  and  5: 


f(j,  ),  however, 
vp 


£  -  E  (a) 


so  that  Equation  14  with  £  given  by  Equation  15  represents  the  desired  consti¬ 
tutive  form  relating  ;  to  e. 
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Exact  Solutions.  For  present  purposes  we  will  consider  solutions  for 
Equation  14  as  if  the  model  represented  a  material  test  specimen  with  either 
a  specified  strain  loading  history  or  stress  loading  history,  and  we  seek 
the  corresponding  response  history. 

To  achieve  an  exact  solution  we  restrict  ourselves  to  linear  functional 

forms  for  k(e  ),  1(f),  and  E(e  )  as  follows: 
vp  e 

k(s  )  *  a  +  H'  e  (If 

vp  y  vp 

<p(f)  -  f/a  ,  (f  >  0)  (1/ 

■  V*  (1! 


with  f  ■  ja!  -  (o  +  H1  e  ) 

11  y  vp 


£  “  J 

vp 


r  ^  | de  j 


where:  E  *  elastic  modulus  (positive  constant). 


j  »  initial  yield  stress  (positive  constant) . 


H'  •  plastic  hardening  modulus  (positive  constant) . 


viscoplastic  strain  hardening  measure. 


Eurtner,  if  viscoplastic  straining  is  monotonic  (i.e.,  £  does  not  change  sign), 

then  e  =  £  1 ,  or  by  Equation  15,  £  *  ■  e  -  o/E„; .  Thus,  the  constitutive 

vp  ■  vp  ■  /  vp  0 

relation  (Equation  14)  becomes  the  following  first  order,  linear  ordinary  dif¬ 
ferential  equation  when  f  >  0: 


6  +  (E  +■  H')o  •  E.(e  +  H'  e  +  r -sgn(loading)  )  (21) 

v  u  u  o 

y  y 

For  the  case  of  stepped  strain  loading  causing  viscoplastic  flow,  we 
specify: 


8 


(22) 


e(t)  -  Be  h(t) 

y 

where  z  ■  a  /E.  ,  initial  yield  strain 

y  y  0 

3  >  1  ,  loading  magnitude 

fo  t  <  o" 

h(t)  =  <  ,  ,  heavyside  function 

I  1  t  >  0 

v_ 

The  solution  to  Equation  21  is: 

j(t.)  =  a  (a  +  (3  -  a)  exp (-y (Eq  +  H')t/a^)) 

where  a  -  (1  +  3H7Eq)/(1  +  H'/Eq) 

consequently,  :  >  a  >  l 


(23) 

(2A) 

(25) 


(26) 

(27) 

(28) 


Some 

(a) 


behaviorial  observations  from  the  solution  are: 

The  largest  stress  occurs  instantly  at  t  =  0+,  a(0)  «  Bc^  ,  which  is 
purely  an  elastic  response,  i.e.,  no  time  for  viscoplastic  flow  to 


occur . 

(b)  The  smallest  stress  occurs  when  t  ■+■  °°,  o(»)  =  ao^,  which  is  the 
elastic-plastic  solution,  i.e.,  steady  state  solution,  *  0. 
(c;  The  rate  of  stress  decay  from  Bo  to  ao  is  controlled  bv  the 

y  y 

fluidity  parameter,  v.  For  large  values  of  v,  o(t)  approaches 
i:j  more  rapidly.  In  the  limit  as  y  +  ”,  o(t)  =  ac  for  all 

y  y 

t  0,  inferring  the  absence  of  the  viscous  element. 

;d)  For  nonhardening,  H'  =  0,  all  the  above  statements  are  valid 
therein  a  =  1. 


The  complimentary  solution  to  Equation  21  for  a  stepped  stress  loading, 
ie fined  by; 


:it) 

is  given  by: 


s(t) 


Bo  h(t) 


B)  exp(-vH't/o  )) 

y 


(29) 


(30) 
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where 


(6  -  1)  +  8 


(31) 


consequently,  a  >  8  >  1 


Some  observations  from  the  solution  are: 

«4» 

(a)  The  smallest  strain  occurs  instantly  at  t  ■  0  ,  e(0)  *  8e  ,  which 
is  purely  an  elastic  response. 

it 

(b)  The  largest  strain  occurs  when  t  ■*  *>,  e(<»)  *  a  t  ,  which  is  the 

y 

steady-state,  elastic-plastic  solution. 

★ 

(c)  The  rate  of  strain  increase  from  Be  to  a  e  occurs  more  rapidly 

y  y 

when  the  fluidity  parameter  y  is  large. 

(d)  When  the  hardening  parameter  becomes  small,  H'  -  0,  we  have 

★ 

a  -►  ®  so  that  e(“)  -*•  «  (as  expected).  In  the  limit,  H'  ■  0, 
Equation  30  becomes,  e(t)  =*  Se^  +  y(3  ~  l)t. 


Numerical  Algorithm.  Beginning  with  the  rate  form  of  Equation  13, 


a  *  E(e  -  e  ) 
vp 

we  integrate  over  one  time  step  from  time  t  to  t  ,  to  get: 

n  n+1  6 


Ao  *  E(Ae-  Ae  ) 
vp 

where,  to  ■  jn+^  -  on 


t  +  At 
•  n 


and,  a  -  a(c  ),  (similarly  for  £  and  e  )  (36) 

n  vp 

Since  is  given  by  Equation  10,  we  will  approximate  Ae^  by  a  Crank-Nicolson 
(12)  time  integration,  i.e., 

Ae  -  At((l-9)en  +  9en+1)  (37) 

vp  vp  vp 


where  6  *  adjustable  integration  parameter  (0  <  9  <  1) 
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A 


If  we  choose  0  =  0,  Equation  37  is  equivalent  to  a  simple  forward  difference 

scheme  (Euler  Method) .  This  is  also  called  an  explicit  method  because  we 

estimate  Ac  based  only  on  information  (e*1  )  at  the  start  of  the  time  step, 
vp  vp 

As  a  consequence.  At  is  restricted  in  size  to  avoid  numerical  instability  (13). 
Alternatively,  if  we  choose  9  >  0,  the  method  is  called  implicit  since  Ae^ 
depends  on  information  at  the  beginning  and  end  of  the  time  step,  generally 
requiring  an  iterative  solution  procedure  within  the  time  step.  For  0  ^  1/2, 
the  implicit  methods  are  unconditionally  stable  (14)  so  that  the  choice  of 
At  is  governed  by  desired  accuracy,  not  stability. 


Returning  to  Equation  34,  with  Ac  replaced  by  Equation  37  and  performing 
_1  VP 

the  E  operation,  we  have  a  numerical  version  of  the  viscoplastic  constitutive 
relation  which  becomes  exact  as  At  0, 

E_1  (Aj)  =  Ac  -  At ( (1  -  0)£n  +  9£n+1)  (38) 

vp  vp 

.Also  we  have  the  side  relations  valid  at  all  time; 


=  -r  5(f)  sgn(  j) 
P 

f  =  '  -  k(c  ) 

vp 


(39) 

(40) 

(41) 


We  now  consider  an  iterative  solution  procedure  for  Equation  38  when  strain 

loading  is  specified,  i.e.  Ac  is  known  at  each  step.  At  time  t  the  quantities 

n  .  n  n 

5  ,  c  and  are  known,  and  che  objective  is  to  determine  these  quantities 

■■■  time  t  For  brevity,  we  assume  the  elastic  modulus  is  constant  at  least 

*  ,  -1  , ,  .  -1  n+1  -1  n 

within  rr.e  time  step  so  that  E  (Ac)  =  E  '  -  E  o  .  Thus,  regrouping 

I'. nation  ’S  with  unknowns  on  the  left  we  have  the  form: 


P(: 


n+1. 


n 

q 


where 


?(an+1) 


e-1 ja+1  +  Ate 


.  n+1 
’"vp 


(42) 

(43) 


n 

q 


Ac 


,  /,  ,\.n  _™1  n 

it(l-9)c  +  E  i 
vt> 


(44) 
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Table  1.  Solution  algorithm  for  one-dimensional  viscoplasticity 
with  strain  loading. 


,  „ .  n  n  n  -n  n  „,n  .  ^  .  ,  n+1 

1.  Given:  e  ,  a  ,  e  ,e  ,  P  ,  P  (at  time  t  )  and  e 

vp  vp  n 


2.  Time  loop:  n  -*•  n+1  up  to  nmax 


/  \  n+1  n  ,  .n  1  n 

(set)  q  3  £  -  e  -  it(l-S)e  +  E  a 

vp 


Iteration  loop:  i  =  1,2,  ...  imax 
(solve)  6oi  *  (P'S  *  (q  -  P*) 
(update)  =  a*  +  Sa^" 


i+1  n+1  -1  i+1 

:  =  -:  -  E  3 


,i+l  r,  i+1  i+1 

f  =  f(a  ,  £  ) 

vp 

.i+1  , ci+l.  ,  i+1, 

-t  =  Yt»(f  )  sgn(o  ) 

i+1  „-l  i+1  ,  .i+1 

P  =  E  o  +  itee 

vp 

i+1  -1 .  3  . . i+1.  * 

’’  =  E  +  it9  —  ( £  ) 

JO  vp 


u.  Repeat  iteration  (step  3)  unless  one  of  the  following  conditions 
are  satisfied: 

(a)  3  =  0  (explicit  solution) 


n  i+1 

(b)  f  and  f  <0  (elastic  space) 


(c)  ,  so  <  tolerance  (converged  solution) 


(d)  i  >  imax  (maximum  iteration  limit) 


5.  Print  results,  return  to  step  2  if  n  <  nmax 


b .  End . 


vp  vp 


_1  u.  ,  i+lx  _ •  i+l. 
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MULTI-DIMENSIONAL  VISCOPLASTICITY 

We  now  consider  a  general  viscoplastic  constitutive  model  for  continuum 
bodies  with  six-dimensional  stress  and  strain  vectors  ordered  as: 

°  "  <ail*  °22’  °33*  al2’  a13*  a23>T 

”  <Cell’  £22*  e33*  y12’  y13’  y23^> 

In  point  of  fact,  the  general  viscoplastic  formulation  is  only  a  small  extention 
of  the  one-dimensional  formulation,  so  that  the  following  development  is  pre¬ 
sented  with  minimal  additional  discussion. 


(49) 

(50) 


Constitutive  Development.  Two  fundamental  assumptions  are  as  follows: 


D(c  ) 

_  .  e 


(51) 


-  e  _vp 


(52) 


where 


€  *  elastic  strain  vector 

_  e 


c  ■  viscoplastic  strain  vector 
_vp 


£  *  elastic  constitutive  matrix  (operator) 

(and  D  .  ,  ,  . 

—  is  assumed  to  exist) 


The  viscoplastic  strain  rate  is  assumed  to  be  of  the  Perzyna  type  (7)  defined  by 
Yd> ( f )  m  (53) 


-VP 


rt( f),  f  >  o 


where. 


>(f)  -t 


viscous  flow  function 


[_°,  f  <  0 


f(o,  e  )  *  yield  function 
-  -vp 


(54) 

(55) 
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?  '(£  held  constant) 
vp 

and,  y  *  fluidity  parameter 


vector  gradient  of  f  w.r.t. a  (56) 


Note  the  above  relationships  are  identical  in  form  to  the  one-dimensional  relation¬ 
ships  wherein  m  replaces  sgn(a) ,  defining  the  vector  direction  of  c  . 

The  yield  function  f  may  be  taken  as  any  valid  plasticity  yield  function 
with  the  following  understanding.  In  six  dimensional  stress  space  all  states 
of  j  giving  f (a,  e  )  =  0  forms  a  "static  yield  surface"  (e.g.  a  classical  plasti¬ 
city  yield  surface),  all  states  of  j  giving  f(a,  e  )  =  3  >  0  forms  a  "dynamic 

yield  surface",  and  all  states  0  such  that  f(j,  £  )  <  0  implies  a  is  in  the 

'  -  -vp 

elastic  domain,  i.e.,  o(f)  =  0  so  that  -;  =  0. 

’  -vp 

Two  typical  forms  of  the  viscoplastic  flow  function,  ■> ,  given  by  Perzyna  are: 

/  r-\* 


:(f)  =  exp  (-f-)  -  1  (59) 

ro 

where  N  =  exponent  (material  parameter)  • 

f^  =  some  normalizing  parameter  so  that  j  is  dimensionless 
(material  parameter) 

equation  is  directly  analogous  to  the  form  adopted  in  the  one-dimensional 
cevelopr.ent  and  is  apparently  the  most  popular  form  for  soils.  The  second  form 
(equation  59)  is  sometimes  used  for  metals.  A  series  generalization  of  these 
two  forms  mav  be  written  as: 


. a.  (i~) 
J=1  J  r0 


f  J 

3  (exp  (— )  -  1) 

J=1  J  "o 


where  A.,  3.  are  additional  material  constants. 
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Having  established  the  fundamental  assumptions,  the  viscoplastic  constitu¬ 
tive  relationship  (i.e.  a  set  of  differential  equations  in  time  relating  o  and 
e)  can  be  readily  derived  by  first  combining  Equations  51  and  52  in  rate  form 


to  get: 


a  «  D(e  -  i  ) 
-vp 


Upon  replacing  e  with  Equation  53,  the  final  constitutive  form  is: 
-vp 


d  »  D(e  _ 


Here  it  is  understood  that  if  f  is  a  hardening  yield  function,  then 

is  an  argument  (i.e.  f  »  f(a,  e  )).  But,  since  e  *  e  -  D-1(a),  Equation  63 

-vp  _vp  .  ” 

implicitly  infers  the  desired  (a,e)  relationship. 

A  numerical  solution  strategy,  paralleling  the  one-dimensional  procedure, 
is  given  next . 

Numerical  Algorithm-  As  in  the  one-dimensional  case,  we  seek  a  solution 
a(c)  when  e(t)  is  specified,  or  vice  versa.  In  either  case,  we  integrate  Equa¬ 
tion  62  over  one  time  step  tn  to  tn+^  and* denote  the  result  by: 


Ao  ■  D(Ae  -Ac  ) 

“  -  .vp 

t  +  At 

,  ,  n+l  n  f  n  .  , 

where  Act  *  a  -  a  *  ;  u  dt 


-  "  - ( c n }  ’  !(tn)*  °r  SUn}  (66) 

Quantities  3°  are  presumed  td  be  known  ac  the  start  of  each  step. 

Using  the  previously  described  Crank-Nicolson  time  integration  scheme  for 

approximating  Ac  ,  we  have  (0  <  6  <  1) : 

-vp  - 


Ae  -  At ( (1  -  9)  £n  +  9£n+1) 
-Vp  .vp  ,vp 
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Inserting  the  above  into  Equation  64  (after  performing  the  £  operation) 
an  incremental  approximation  for  the  viscoplastic  constitutive  relationship 
is  obtained  (becoming  exact  as  At  0); 

J2-1(Ao)  =  Ae  -  At ( ( 1  -  9)en  +  9en+1)  (68) 

-  -  _vp  _vp 

Implied  in  the  above  are  the  viscoplastic  model  assumptions  valid  for 
all  time: 

£  =  Yb(f)m  (69) 

~vp 

-1 

tvp  =  f  -  °  <70> 

For  the  case  when  strain  loading  is  specified.  Equation  68  is  regrouped  with 

unknowns  at  time  t  ,  on  left; 

n+1 

?.:rn+1)  =  qn  (71) 


n+1. 

where  ?(j  ) 


-l,n+l  .  .n+1 

D  -  +  S  t  3  i 


qa  =  1c  -  At (1  -  9) -:n  +  D_1on  (73) 

-  -  -vp  -  - 

Here  it  is  assumed  the  elastic  matrix  is  constant  (at  least  within  the  time  step) 
If  implicit  integration  is  employed  (9  >  0) ,  Equation  71  represents  a  set 


c-  nonlinear  algebraic  equations  for 


Using  a  Newton-Raphson  procedure  we 


n+-  i  .  i  i  .  '  .  .  n+1  ,  .  i. 

..ec  :  =  o  +  to  wnere  z  is  some  estimate  or  z  and  to  is  a  correction. 

.'-t  approximation  for  5a1'  is  achieved  by  expanding  P(o1  +  So1)  in  a  first  order 

Taylor  series,  and  solving: 


_  ,  ,  _  i,  .  i  n  „  ,  i. 
L  ■  '  >  'O  »  q  -  P  (z  ) 


wr.er:  the  Jacobian  matrix  is: 


2.'  (:*) 


SPfo1) 


and  ?(oS  =■  2  +  At9s^ 

-  -  -vp 
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The  process  is  repeated  with  ■  a*  +  5q^  to  get  5ai+^  and  so  on,  until 

1 5<?  |  -  0.  Table  2  summarizes  the  algorithm  for  strain  loading. 

For  the  case  of  stress  loading,  a  similar  solution  procedure  is  readily 
derived  to  iteratively  determine  gn+^  *  e*  +  <Se*,  i.e.; 


Z'  (e1)6ei  -  qn  -  P(e  ) 

where  q11  -  D  ^Ao  +  At(l-6)en 

,  "  -Vp 

PCe1)  -  e1  -  0ACE1 
-vp 

ruh  -  3?(cl) 

3e 


(77) 

(78) 

(79) 

(80) 


The  solution  procedure  for  stress  loading  is  summarized  in  Table  3.  Note  the 
solution  procedure  for  stress  loading  is  very  similar  to  the  strain  loading 
procedure  wherein  q,  P,  and  are  replaced  by  q,  P,  and  . 


The  complete  Jacobian  matrix  Z'  (or  for  stress  loading)  is  generally 

nonsymmetric  when  hardening  is  included  in  the  yield  function.  This  may  be 

avoided  by  not  including  contributions  to  the  Jacobian  matrix  arising  from 

partial  derivatives  of  f(o,  e  )  w.r.t.  e  as  illustrated  below. 

•  -vp  -vp 

Equation  75  defines  the  complete  Jacobian  with  the  understanding  that 
the  hardening  argument  in  f(o,  ■  )  is  to  be  replaced  by  e^  «*  g  -  D  ^o,  so 

that  tne  complete  Jacobian  matrix  is: 


h 


9m 


(symmetric) 


(81) 

(82) 

(83) 

(84) 

(85) 
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Table  2.  Solution  algorithm  for  general  viscoplastic  model  with 
strain  loading. 


n+1 


1. 

_ .  n 

Given :  e 

n  n 
,a  ,e  , 
’-vp 

-n  pn  p  > n  r . 

-vp’  ?  ’  1  ° 

2. 

Time  loop: 

n  -*■  n+l 

,  up  to  nmax 

.  .  n 

(set)  q 

,  n+l 

■  \  § 

-  en)  -  At ( 1  - 

3. 

Iteration 

loop :  i 

=  1,2,  ...  imax 

(solve)  : 

So1 

n  T,  i 

=  q  -  P 

(update)  : 

i+1 

a 

i  ,  .  i 
a  +  oq 

i+l 

P  33 

-vp 

n+l  -1  i+l 

c  "Si 

£  i+l 
*.  = 

i+l  i+l. 

'<?_  •  > 

_  i+l 

;  3fi+1 : 

J  J 

: i+1 
':vp 

,  .i+l"  i+l 

■' :  ( r  )  m 

?i+1  = 

At3c 

.n 

P 


-1  n 


vp 


3Fi+1 

3o 


Repeat  iteration  (step  3)  unless  one  of  the  following  is  satisfied 


(a)  3=0,  (explicit  integration) 

(b)  fn  and  <  0  ,  (elastic  space) 

(c)  ij1'  tolerance,  (converged  solution) 


(d)  i  '  imax  ,  (iteration  limit) 
Print  results,  return  to  step  2  if  n  <  nmax 
End 


Symmetric  form  of  Jacobian  matrix  is  (see  Equation  87) 


-1  T 

D  +  ylAt [o’  m  m‘  +  3  m’ i 
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Table  3.  Solution  algorithm  for  general  viscoplastic  model  with 
stress  loading. 


1. 

nn.nn_nn.  .  .  ,  n+1 

Given:  e  ,  a  ,  g  ,  £  ,  p  ,  P  (at  time  t  ),  and  a 

vp  -vp  “  n 

2. 

Time  loop 

n  -*■  n+1 

,  up  to  nmax 

(set)  q 

f  n+ 1  n.  ,,  v.n  n 

*  0  (a  -  a  )  +  At  (1  -6)  e  +  e 

-vp 

3. 

Iteration  loop:  i  = 

1,2  ...  imax 

(solve) : 

~ ,  i  .  i 

P/  6e  = 

~n  li 
q  -  ? 

(update) : 

i+1 

£  * 

i  ,  .  i 
£  +  <5e 

i+1 

e  * 

-vp 

i+1  _-l  n+1 

£  -  £  a 

i 

1 

fi+1  = 

i+1  _  , 

f(0n+1,  e^1) 

r  ~  -  ^ 

3fi+1  N 

1 

3o 

1 

.i+1 

e  ■ 

-vp 

s_./.ri+lv  i+1 
y<S(f  )  m 

1 

Pi+1  , 
r  * 

ti+1  -  At9ei+1  ! 

’  t-  Vvp  ! 

1  3Pi+1] 

1  -  1  *  ! 

:  a?  i  ! 

~  ~ *  i 

4. 

» 

Repeat  iteration  (step  3)  unless  one  of  the  following  is  satisfied:  j 

(a)  6=0,  (explicit  integration) 

(b)  fn  and  fi+^  <  0  ,  (elastic  space) 

/  _  \  -  i+1  :  i _  /_ 


(c)  | 5c  |  <  tolerance  ,  (convergence  solution) 

(d)  i  >  imax  ,  (iteration  limit) 

Print  results,  return  to  step  2  if  n  <  nmax. 

End 


Symmetric  form  of  Jacobian  is  identify  matrix,  P'  =  thus  the  "solve" 
step  is  trivial  and  P'  need  not  be  formed. 
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i 


1 


(non  symmetric) 


(86) 


jm 


‘3c 


-vp 


and  all  the  above  partials  are  taken  with  the  other  variables  held  constant. 

If  f  ■  f(a),  implying  no  hardening,  then  £'  becomes  a  symmetric  matrix 

because  _h  and  are  2ero.  Otherwise  £'  has  nonsymmetric  contributions  from 
T  T 

the  matrices  m  h  and  .  By  retaining  only  the  symmetric  terms  to  define 
we  have: 


-1  T 

£'  =  S  +  Y0At(b'  mm*  +  <p  a' ] 


(87) 


In  a  similar  development  for  stress  loading,  we  find  is  composed  of  the 
identity  matrix  plus  nonsymmetric  matrices  related  to  hardening.  Retaining  only 
the  symmetric  identity  matrix  we  have: 

l'  -  I  (88) 

In  closing  this  section  it  is  emphasized  that  the  symmetric  forms  of  the 
Jacobian  matrices  (Equations  87  and  88)  are  complete  and  correct  for  non  hardening 
yield  functions.  For  the  more  general  case  of  hardening,  the  symmetric  forms 
may  be  used  with  the  consequence  that  the  N'ewton-Raphson  rate  of  convergence 
anc  the  domain  of  convergence  (i.e.  step  size)  may  be  reduced.  Of  coursje  once 
a  converged  solution  is  obcained,  it  is  correct. 

As  illustrated  in  the  last  section  of  this  report,  the  inverse  of  the 
Jacobian  matrix  becomes  the  "constitutive  matrix"  in  forming  finite  element 
stiffness  matrices.  Thus  the  motivation  for  retaining  the  symmetric  form  is 
apparent  with  regard  to  equation  solving. 
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SPECIALIZATION  OF  VISCOPLASTIC  MODEL  TO  CAP 7 5 

Functional  Forms.  The  foregoing  viscoplastic  algorithm  presented  in 

Table  2  (or  Table  3  for  stress  loading)  is  specialized  by  identifying  the 

functional  forms  for  elastic,  viscous,  and  plastic  components,  i.e.  $(f), 

and  f(o,  e  )  along  with  the  hardening  function  associated  with  f. 

-  -vp 

In  program  VPDRVR  (Appendix)  the  following  forms  are  included: 

(1)  An  isotropic  elastic  matrix  is  formulated  with  a  bulk  and  shear 
modulus  (K  and  G)  whose  values  may  change  between  load  steps 
but  are  assumed  constant  within  the  load  step,  i.e.; 

£  -  £(K,G)n  (89) 

(2)  Two  options  for  4> ( f )  are  (see  Equations  58,  59)  : 

r  (f/fo>N 

$(f)  *  ^  (or)  (90) 

"  exp  (f/fQ)N  -  1 

therefore, 

,Vo(f)/f0 

*'(f)  -  <  (or)  # 

L  N(j(f)  +  1)  fN_1/f0N  (91) 

(3)  A  somewhat  general  form  for  f(o,  s  )  applicable  to  CAP75  as  well 
as  most  other  plasticity  models  is: 

f(o,  £  )  «  g  (J  ,  e  )  +  g  (J  )  (92) 

-  -vp  1  1  -vp  °2  2 

where  g,  =  specified  function  of  J,  and  t  where  e 

1  1  -vp  -vp 

implies  hardening. 

g?  *  specified  function  of  J?. 

•  first  invariant  of  stress. 

J,  =  second  invariant  of  deviator  stress. 
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In  a  later  section  specific  forms  for  g^  and  are  given  for  the  CAP75  plasticity 
model  where  we  allow  and  g^  Co  be  specified  differently  in  different  regions 
of  J^,  J7  space  (e.g.  failure  surface  and  cap). 

The  advantage  of  the  above  form  for  f  is  that  the  gradient  vector  m  and 
the  Jacobian  matrix  £/  (symmetric  form)  can  be  established  once  and  for  all 
in  terms  of  g^  and  g^  and  their  derivatives  as  shown  below. 

Beginning  with  the  stress  invariant  definitions: 


all  +  °22  +  a33 


(93) 


J2  ‘  1  <SU  + 

2  2  2  2 

S  +  S  +  2o  +  2a 

22  33  12  13 

+  24> 

(94) 

where  S  .  m 

n 

aii  ~  V3’  1  =  1*2’3 

(no  sum) 

(95) 

we  have  for  a  =  -If /3c 

: 

?  =  b  +  §; 

3 

(96) 

where  = 

(97) 

’*2 

>  ^  7 

(98) 

b  = 

=-'111000 

T 

•> 

(99) 

a  = 

!V*  ■  *  Su  s22  S33 

T 

2(7  2(7  > 

12  13  23 

(100) 

Anc  after  some  manipulation,  2.'  (Equation  87) 

may  be  written  as: 

T 

L  =  2'"  + 

t  (X^  [a  a  *  ]  -t-  [baT  + 

a  bT]  +  X  [b  bT]  +  X, 

~  3  -  ~  4 

la 

[t-D 

70 

(101) 

where  = 

(102) 

X2  = 

®  ®1^2 

(103) 

23 


(104) 


X3  ”  9'  ^8p2  +  ^81 

X  -  <dg'  (105) 

4  2 


(106) 

(107) 

(108) 


Thus  to  change  yield  functions,  we  merely  specify  new  functions  for  g^  and 
with  g^,  g ”,  g^,  and  g^ •  Of  course,  hardening  parameters  must  be  maintained 
separately  in  updating  the  g,  function. 

CAP75  Yield  Function.  We  now  specialize  Equation  92  for  a  particular 
version  of  the  CAP75  plasticity  soil  model  (1) .  Figure  2  shows  the  CAP75 
yield  function  which  is  composed  of  a  nonhardening  "failure  surface"  f  and 
a  hardening  "Cap  surface"  f  .  The  failure  surface  is  defined  as: 


(1C9) 

(110) 

(111) 


yote  the  so-called  "von  Mises  transition"  is  not  included  in  this  development, 
and  the  ad  hoc  "tension  cutoff"  technique  (1)  is  not  suitable  in  a  viscoplastic 
context  (see  Summary) . 
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The  cap  surface,  forming  an  ellipse  quadrant,  is  defined  as: 


fC  “  8C1^J1»  e) 

+  *c2(V 

(112) 

where,  gc 

(-(X  -  L)2  +  <JX  -  U2)/£0R2 

(113) 

8Co 

Vf0 

(114) 

and,  e  *  volumetric  viscoplastic  strain  hardening  function 

R  *  ratio  of  principal  ellipse  radii 

L  *  current  location  on  axis  of  cap-failure  intersection 
X  »  current  location  of  cap  intersection  on  axis 
(X<  L) 

f  •  normalizing  constant  with  units  of  stress. 

Note,  f  is  a  "squared”  form  of  the  cap  surface  presented  by  Sandler  (1) . 

This  is  necessary  because  during  viscoplastic  flow  (J^<  X)  we  must  have 
fc  >  0  which  is  not  possible  with  the  Sandler  form.  Accordingly,  f^  is  used 
to  retain  the  units  of  stress  for  f  . 

Continuity  of  the  hardening  cap  and  failure  surfaces  (static)  are  main¬ 
tained  fcy  the  relationship; 

X  -  L  +  R  g  (L)  (115) 

F1 

which  establishes  the  current  ellipse  quadrant  of  the  cap  surface.  Typically, 
the  ellipse  ratio  R  is  assumed  constant  but  may  be  a  function  of  L. 

CAP 7 5  Hardening  Function.  Hardening  of  the  cap  is  given  by  a  relationship 
between  X  and  e; 

e  -  W(eDX  -  1)  (116) 

where  W  and  D  are  positive  material  constants. 

The  strain  hardening  function  e  is  a  measure  of  accumulated  volumetric 
viscoplastic  strain.  For  both  rocks  and  soils,  e  can  only  increase  negatively 
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during  viscoplasCic  loading  of  Che  cap,  inferring  X  and  L  grow  negatively 
expanding  the  cap  during  hardening.  Thus  the  incremental  update  of  e,  when 
Jx  a,  is: 


n+1 


e  n  +  min(0,  Aw) 


where 


,  n+1  n 

Aw  ■  w  -  w 


vp 


+  e 


11 


vp 


+  £ 


22 


vp 


33 


(117) 

(118) 
(119) 


For  soils  only,  e  is  allowed  to  be  updated  with  positive  increments  of 
Aw  during  viscoplastic  loading  of  the  failure  surface,  inferring  that  X  and  L 
become  less  negative.  This  retracts  the  cap,  however  the  retraction  is  limited 
such  that  L  <  J^.  Sandler  claims  this  is  tantamount  to  kinematic  hardening, 
not  strain  softening.  The  purpose  of  this  procedure  is  to  limit  excessive 
dilatancy  for  soil  models  (i.e.,  subsequent  compression  loading  will  activate 
the  cap  sooner).  Thus  the  incremental  update  of  £  for  soils,  when  >  L,  is: 


£  »  min  (£a,  c^)  (120) 

where  £  ■  e11  +  max  (0,  Aw)  (121) 

a 

c  -  W(eDXc  -  1)  (122) 

C 

Xc  -  Jl  +  Rgp  (Jj)  (123) 


Irs  the  abo/e,  insures  that  L  <  J.  .  Also,  e  is  never  allowed  to  become 
c  —  1 

greater  roan  its  initial  value. 


Table  A  summarizes  the  implementation  of  the  CAP75  plasticity  model.  The 
procedure  in  Table  A  should  be  viewed  as  a  subroutine  which  internally  main- 
'.air.:  the  hardening  parameters  e,  X,  and  L.  The  subroutine  outputs  g^,  g| ,  g^ 
and  g7>  &2>  §2  which  allows  calculation  of  f,  m,  and  J»*  by  Equations  92,  96, 
and  101,  respectively.  These,  in  turn  provide  the  necessary  ingredients  in 
Table  2  (or  Table  3  without  £' )  for  the  general  viscoplastic  algorithm. 
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Input  initialization:  Initial  location  of  X  specified  as  X 

check:  X°  <  g_  (0) 

F1 

set:  *  W(exp(DX^)  -  1) 

solve:  *  X^  -  Rg  (L^) 

F1 


For  each  iteration  of  each  load  step:  Update  e,  X,  and  L  and  compute 
81’  81’  81  and  8 / ’  82*  8 2  based  on  current  values  of  J^,  J2  and  Aw. 


J1  1  L  ? 


Failure  Surf .  „  / 


For  soil  onxy : 

f  +  max(0,  Aw) 
e  =  min  l  W(exp(DX  )  -  1) 

L  e° 


-  J  +  RgF  (J .) 


Yes 

\  Cap  Surf. 


t  *  en  +  min(0.  Aw)  j 


ln(e/W  +  1)/D 


X  -  RgF  (L) 
F1 


ail  Surf. 


EXAMPLE  RESULTS  AND  DISCUSSION 

In  this  section  the  CAP7S  viscoplastic  model  is  studied  for  tvo  loading 
conditions;  uniaxial  strain  and  triaxial  stress.  In  both  cases  the  plasticity 
parameters  are  taken  from  Sandler  (1)  representing  McCormick  Ranch  Sand. 

Table  S  gives  the  values  of  the  model  parameters  used  for  both  loading  con¬ 
ditions  along  with  three  values  for  the  fluidity  parameter. 

The  intent  of  these  examples  is  to  verify  the  viscoplastic  algorithm  and 
to  illustrate  the  effect  of  the  fluidity  parameter,  y,  on  the  time  dependent 
responses.  The  three  values  of  the  fluidity  parameter  (y  *  0.001,  0.01,  and 
0.1)  were  chosen  to  produce  responses  with  varying  amounts  of  viscoplastic 
flow. 

With  regard  to  numerical  time  integration,  each  example  problem  is  solved 
within  It  relative  accuracy  for  five  choices  of  the  Crank-Nicolson  9  parameter; 

9  »  0.0,  0.25,  0.50,  0.75,  and  1.0.  That  is,  for  each  choice  of  9,  the  time 
step  size.  At,  is  repeatedly  halved  until  successive  solutions  agree  within  It, 
thus  providing  an  efficiency  comparison. 

Uniaxial  strain.  Figure  3  shows  the  uniaxial  strain  loading  history  for 
the  strain  component  s^,  all  other  strain  components  are  zero.  Here,  in¬ 
creases  in  compression  at  a  constant  rate,  held  constant,  unloaded  at  a  con- 
:tant  rate,  and  finally  held  constant  so  that  eventually  stress  responses  would 
approach  steady  state  for  all  y.  Note  the  units  of  time  are  not  and  need  not 
be  explicitly  stated  since  they  are  the  reciprocal  of  whatever  units  are  assumed 
for  v . 

Figure  4  shows  the  corresponding  stress  response  for  the  three  magni¬ 
tudes  of  y.  For  the  relatively  large  magnitude,  y  *=  0.1,  the  response  is  nearly 
inviscio  throughout  the  entire  loading  history.  In  other  words,  the  viscous 
element  (see  Figure  1)  has  very  little  effect  in  retarding  the  plastic  response 
cQ  that  the  solution  is  nearly  elasto-plastic.  Indeed,  the  steady  state 
values  agree  exactly  with  the  inviscid  CAP75  plasticity  solution  (17),  thereby, 
providing  a  verification  of  the  viscoplastic  solution.  If  higher  values  of  y 
were  used,  the  response  would  not  change  significantly,  and  the  small  peaks 
at  time  *  1.0  and'  5.5  would  be  flattened,  resulting  in  a  perfectly  inviscid 
so  lution. 
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Table  5.  Viscoplastic  model  parameters 


Elastic  Parameters 

Bulk  Modulus  -  66.67  ksi 

Shear  Modulus  ”  40.0  ksi 

CAP75  Plastic* 

Failure  Surface  (Equation  109) : 

(McCormick  Ranch  Sand) 

i 

i 

i 

i 

1 

\ 

A  -  0.25,  B  -  0.67,  C  -  0.18 

Cap  Surface  (Equation  112) : 

R  -  2.5,  f  -  0.25,  X°  -  -0.1888 

0 

Cap  Hardening  (Equation  116) : 

W  -  0.066,  D  -  0.67 

\ 

Viscous  Parameters 

I 

! 

• 

! 

i 

* 

Flow  Function  (Equation  90a) 

N  »  1,  f  =  0.25 
° 

Fluidity  parameter 

Y  -  0.001,  0.01,  and  0.1 

*  o  -1-1 

Units:  A,  C,  fQ  and  X  ■  ksi;  B  and  D  *  ksi  ,  y  »  time  ; 


R,  W,  and  N  are  dimensionless. 

I 

! 
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Figure  A.  Axial  stress  response 


For  the  relatively  small  magnitude,  y  *  0.001,  Figure  4  shows  a  pronounced 
viscous  effect  such  that  steady  state  conditions  are  not  reached  within  the 
constant  loading  duration,  and  peak  stress  magnitudes  far  exceed  the  maximum 
stresses  achievable  from  the  inviscid  plasticity  solution.  The  rate  of  stress 
relaxation  decreases  with  decreasing  values  of  y.  In  the  limit  as  y  +  0,  the 
response  becomes  perfectly  elastic  in  which  case  the  shape  of  the  o  response 
history  would  be  similar  to  the  loading  history. 

For  reference.  Figure  4  also  indicates  which  CAP75  regions  (i.e.  failure 
surface,  cap  surface,  or  elastic  domain)  are  being  activated  during  the  loading 
schedule.  Lateral  stress  responses  (o ^  ”  ^33)  are  shown  in  Figure  5. 

The  foregoing  solutions  were  independently  obtained  for  five  choices  of  9, 
each  requiring  a  certain  minimum  value  of  At  to  achieve  relative  accuracy  within 
17..  Table  6a  shows  the  required  values  of  At  to  maintain  accuracy  (Note, 
although  not  necessary.  At  was  taken  uniformly  throughout  the  loading  schedule). 
In  general  it  is  observed  that  9-0.5  permits  the  largest  time  step,  in  some 
cases  up  to  four  times  greater  than  9  -  0.0  or  9  -  1.0.  Surprisingly,  the 
influence  of  y  on  At  is  small.  That  is,  within  the  range  0.001  <_  v  _<  0.1,  At  for 
accuracy  changes  at  most  by  a  factor  of  two.  This  is  not  in  conformance  with 
initial  expectations  wherein  it  was  incorrectly  anticipated  that  At  for  accuracy 
would  be  inversely  proportional  to  y.  In  other  words,  it  wls  presumed  the 
dimensionless  product  yAt  would  be  an  accuracy  measure. 

It  is  now  conjectured  that  the  accuracy  requirement  for  At  is  predominantly 
controlled  by  the  loading  schedule,  i.e..  At  must  be  sufficiently  small  to 
adequately  capture  the  abrupt  changes  in  loading  rates.  This  contention  is 
supported  by  the  observation  that  the  greatest  accuracy  deviations  throughout 
the  loading  schedule  occurred  at  changes  in  loading  rates,  primarily  at  time  = 

■*  r\ 

A..  J  • 

Lastly,  the  At  values  in  Table  6a  for  the  conditionally  stable  algorithms, 

-•  *  0.0  and  0.25,  are  controlled  strictly  by  accuracy  not  stability.  Unstable 
values  of  At  were  typically  found  to  be  an  order  of  magnitude  or  more  than  the 
accuracy  controlled  values  at  At.  Unstable  solutions  are  readily  distinguished 

ov  wild  oscillations  of  the  responses. 

Triaxial  Stress.  Figure  6  shows  the  triaxial  stress  loading  schedule. 
Initially,  hydrostatic  loading  (a^  -  j22  *=  j^)  is  applied  at  a  constant  rate 
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Table  6a.  Uniaxial  strain:  required  At  for  accuracy 


Integration 

Fluidity  Parameter 

Parameter 

0 

IIIPB88SBM 

'  ■ 

Y  -  0.1 

0.00 

0.00625 

0.0125 

0.00625 

0.25 

0.00625 

0.025 

0.0125 

0.50 

0.025 

0.025 

0.0125 

0.74 

0.00625 

0.0125 

0.0125 

1.00 

0.00625 

0.0125 

0.00625 

Table  6b.  Triaxlal  stress:  required  At  for  accuracy 


1 

1 

Integration 

Parameter 

; 6 

Fluidity  Parameter 

- ■■■-— - -  1 

Y  -  0.001 

Y  -  0.01 

Y  -  0.1 

i 

0.00 

0.050 

0.0125 

| 

0.00625 

0.25 

0.10 

0.025 

0.0125 

0.50 

• 

0.10 

0.05 

0.025  j 

0.75 

0.05 

0.025 

i 

0.0125 

1.00 

0.025 

0.0125 

0.00625  j 

and  then  held  constant.  Next,  while  holding  the  lateral  stresses  constant, 
the  axial  stress,  o^,  is  increased  at  a  constant  rate,  held  constant,  unloaded 
at  a  constant  rate,  and  finally  held  constant  at  hydrostatic  pressure. 

The  corresponding  axial  strain  response  is  shown  in  Figure  7  for  the  three 
magnitudes  of  y.  As  in  the  previous  example,  the  relatively  high  magnitude, 

Y  *=  0.1,  produces  a  nearly  inviscid  elasto-plastic  response  history.  Accordingly, 
large  amounts  of  plastic  straining  is  observed.  Conversely  for  the  relatively 
low  magnitude,  y  »  0.001,  rapid  plastic  straining  is  impeded  resulting  in  smaller 
strain  magnitudes.  In  the  limit  as  y  -*■  0,  the  strain  response  would  be  purely 
elastic.  Figure  8  shows  the  corresponding  response  histories  for  the  lateral 
strains  *  ^33)  * 

Table  6b  lists  the  minimum  values  of  At  to  achieve  accuracy  within  12  for 
the  five  choices  of  8.  As  before,  9  »  0.5  is  the  most  efficient  choice  for 
maintaining  accuracy  with  maximum  At.  Further  it  can  be  observed  that  the  in¬ 
fluence  of  v  is  fairly  uniform.  That  is,  for  any  value  of  9,  At  generally  de¬ 
creases  by  one  half  for  each  magnitude  increase  in  y. 

Summarizing  both  examples,  the  following  observations  are  noted: 

1.  For  y  >>  0.1,  the  responses  are  inviscid  (elasto-plastic).  For 
Y  <<  0.001  the  responses  are  elastic. 

2.  For  all  y,  9  =  0.5  is  most  efficient  with  regard  to  maintaining 
accuracy  with  largest  At. 

3.  In  the  range,  0.001  <_  y  _<  0.1,  the  fluidity  parameter  has  less 
influence  on  the  required  At  for  accuracy  than  anticipated. 

Loading  rates  and  changes  in  loading  rates  are  presumed  to  be 
the  major  factor  for  the  <-t  accuracy  requirement. 

4.  For  9  <  0.5  (conditional  stability),  the  At  required  for  accuracy 
is  significantly  less  than  the  At  limit  for  stability. 

It  is  emphasized  that  the  above  observations  pertain  only  to  the  example 
problems  presented  herein. 
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Figure 


EXTENSION  TO  FEM  APPLICATIONS 

For  the  sake  of  completeness,  this  section  illustrates  how  to  incorporate 
the  general  viscoplastic  constitutive  model  into  a  quasi-static,  finite  element 
displacement  formulation. 


Using  standard  FEM  notatioo,  the  global  equilibrium  equations  to  be  satisfied 
at  any  time  t  ^  are: 


r  RT  <rn+1  AM 

ly  I  ?  dV 


^+l 


(124) 


with  £  *  JL  u 


(125) 


where 


£  =  strain-displacement  matrix 

-  ■  nodal  point  displacements 

F  =  nodal  force  vector  (body,  traction,  and  point  loads) 
V  »  volume  of  body 


Since  an+^"  is  the  unknown  in  Equation  124,  we  follow  a  procedure  paralleling 
Table  2  far  strain  loading,  the  only  difference  being  that  strain  is  not  a  priori 
specified  but  determined  iteratively  via  displacements.  To^  this  end  we  let; 


n+1 


5c1 


n+1 


+  ae1 


i  i  .  n+1  n+1 

wnere  a  ,  e  =  estimates  or  a  ,  e  at  iteration  i. 


.1,1  .  .  i  l 

,  le  *  correction  to  estimate  o  ,  £  . 


rhus,  the  equilibrium  equation  may  be  written  as: 


3T  5a1  dV  -  6  F1 


with  i F^  »  Fn+^  -  :'v  a1  dV 


Note  when  o*  *  an+\  we  have  5F1  *  0. 


(126) 

(127) 


(128) 

(129) 
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From  Equation  74  we  have  a  relationship  for  dq  : 

P’  6a1  -  (de*  +  e1  -  z*)  -  D-1^1  -  o'1)  -At((l-9)eQ  +  Qi  1  )  (130) 

'  -  vp  ~vp 

Or  more  simply  (noting  ■  e  -  £  ^a)  : 


do*  -  fi1  (Se1  +  dqS 


(131) 


i  .  i  n  .  «  •  n  _  i  . 

where  dq  ■  (e  -  e  )  -At  ((1-9)  e  +  9e  ) 
3  -vp  -vp  -vp  -vp 

m  [p> (a1)]  ^  Jacobian  inverse 


(132) 

(133) 


In  Equation  132,  both  groups  of  terms  defining  6q  are  estimates  of  the 

increment  Ae  when  these  estimates  agree  (converged  solution) ,  we  have  dq*  ■  0. 
-vp 

Upon  inserting  Equation  131  into  Equation  128  with  de1  -  gdu1,  the  iterative- 
global-equilibrium  equation  becomes: 


K1  du1  -  dF1  -  5Q1 


(134) 


where  *  /v  BT  g1  £  dV 


dQ1  •  g1  gL  dq1  dV 


(135) 

(136) 


,  ,  -  . „  .  , .  x  ,  i  n+1  i  n+1 

As  tne  number  or  iterations  increases  (i  -*■  ®) ,  we  have  a  =  a  ,  e  “  z  , 

=  cn+\  etc.,  as  well  as,  dF*  3  dq1  3  dQ^  3  du1  *  0,  (i.e.  a  converged 
■vp  -vp  .  3 

solution) . 


Table  7  summarizes  the  numerical  FEM  algorithm  for  a  time  step  from  time  t^ 

to  t  .  First  iteration  estimates  (i  *  1)  are  based  on  previously  known  values 

n+  ,  1  n  1  nln  1  n 

at  time  t  so  that  we  have  a  •  a  ,  e  •  t  ,  £  *  £  ,  *  &  ,  etc.  Accordingly 

from  Equations  129  and  132,  we  have  5F^  3  Fn+^  _  F°,  and  dq^*-AteJ^. 

Several  variations  of  the  solution  strategy  in  Table  7  are  easily  established. 
First,  for  the  explicit  method  (9  »  0)  we  have  g1  »  D  resulting  in  linear  equations 
chat  convergence  in  one  iteration  but  at  the  expense  of  reduced  At  to  control 
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Given  at  Gauss  Points:  on,  £n,  £  n  ,  en  ,  and  £n  (or  reconstruct), 
n  -vp  -vp 

along  with  global  IL  • 

Time  loop:  n  ■+  n+1,  up  to  nmax. 

(set:)  6F1  -  F(t  j.)  ~  F(cn> 

6qX  =  -At i  n 
-vp 

5Q1  -  /v  S=T  £n  Sq1  dV 

Iteration  loop:  i  =  1,  2,  ...  imax. 

(solve)  :  J^du*  ■  6F  -  dQ'*' 

.  i  _T  .  i 

de  =  B  du 

,  .  ,  i+1  i  .  .  i 

(update):  e  »  e  +  oe 


j1+1  =  a1  +  ^(de1  +  dq1) 

i+1  i+1  -1  i+1 

e  =  e  -Da 


.  i+1  . /<ri+l.  i+1 

£  *  y<K  f  )  m 

,vp 


i+1 

£ 

ir1*1]-1  * 

5Fi+1  = 

,  .  .  T  i+1  i+1  ... 

r(t„+i>  ■',*£  z  dv 

dqi+L  - 

s1+1  -  En  -  AcCCl  - 9)e  n  +  9e 

.vp  .vp  .vp 

5Qi+1  - 

/v  gj  £1+1  5q1+1  dV 

£i+1  - 

/v  BX  £1+1  a  dV 

Repeat  iteration  (step  3)  until  convergence  (see  Table  2) . 
Print  results,  return  to  (step  2)  if  n  <  nmax. 

End. 


-1  -  T 

D  +  YSAtf'j'm  m  +  Pja'  ] 


(see  Equation  87) 
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accuracy/stability.  A  second  variation  is  with  regard  to  modified  or  quasi 
Newton-Raphson  methods  (15) .  Here  is  held  constant  (modified  Newton-Raphson) 
or  approximately  updated  in  its  triagularized  form  (quasi  Newton-Raphson)  for 
the  purpose  of  reducing  solution  time  per  iteration  but  at  the  expense  of  slower 
convergence  rates.  Lastly,  if  the  iteration  process  is  terminated  prior  to 
convergence  (or  if  the  one-step  explicit  algorithm  is  used)  the  out-of-balance 
force  <5F^+^  may  be  added  to  the  applied  force  F  in  the  next  time  step  as  a 
"load  correction"  to  reduce  error  accumulation. 

For  dynamic  problems  the  acceleration  vector  can  be  numerically  integrated 
by  any  suitable  scheme  (e.g.  Newmark  Beta-method)  independently  of  the  Crank- 
Nicolson  scheme  used  for  the  viscoplastic  constitutive  model.  In  such  cases, 
it  may  be  presumed  that  At  controlling  accuracy  would  not  be  smaller  than  At 
required  for  same  problem  with  inviscid  plasticity  (i.e.  CAP 7 5  without  viscous 
effects).  This  is  because  the  viscoplastic  response  at  any  instant  is  bracketed 
oetween  a  purely  elastic  response  and  an  inviscid  plastic  response.  Thus,  it  is 
reasonable  to  assume  that  At  for  a  dynamic  viscoplastic  solution  need  not  be  less 
than  At  for  a  dynamic  inviscid-plastic  solution. 
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SUMMARY  AND  RECOMMENDATIONS 


The  viscoplastic  formulation  and  numerical  algorithm  developed  herein 
provides  a  general  format  for  incorporating  various  plasticity  models  into  a 
Perzyna-type  viscoplastic  constitutive  relationship  suitable  for  FEM  appli¬ 
cations.  In  particular,  the  CAP75  viscoplastic  model  illustrated  in  this 
report  and  contained  in  the  program  VPDRVR  (Appendix)  appears  to  have  a 
sufficient  generality  to  faithfully  represent  the  time-dependent  behavior 
of  many  geological  materials  over  a  vide  range  of  loadings. 

The  VPDRVR  program  has  been  extensively  checked  and  verified  including 
cross-checks  with  the  CAPDRVR  program  for  steady  state  elasto-plastic  solutions, 
as  well  as,  self-checking  time-dependent  inverse  solutions.  Inverse  solutions 
are  obtained  by  taking  the  stress  responses  from  a  strain  loading  problem  and 
using  them  for  a  stress  loading  problem  whose  strain  response  should  match  the 
original  strain  input  (or  vice  versa) .  These  severe  self-checking  tests  demon¬ 
strated  that  the  algorithm  is  accurate  and  working  admirably.  Furthermore, 
che  architecture  of  the  VPDRVR  program  permits  relatively  easy  addition  of  new 
elastic,  plastic  and/or  viscous  functional  forms. 

Future  research  is  needed  to  bring  this  work  to  full  fruition.  Two  major 
areas  are  (1)  experimental  verification  and  parameter  identification,  and 
(2)  numerical  studies  for  optimizing  ef ficiency/accuracy  of  the  FEM  algorithm. 

The  former  area  should  be  addressed  first  to  establish  the  capabilities 
and  limitations  of  the  viscoplastic  model.  Existing  experimental  data  for 
time-dependent  behavior  of  soils  and  rocks  (2, 3, 4, 5)  is  primarily  for  slow 
rates  of  loading.  Hence,  it  is  strongly  recommended  to  obtain  additional 
experimental  data  for  rapid  loading  rates.  Concurrent  with  experimental  veri¬ 
fication,  is  che  need  for  developing  feasible  parameter  identification  tech¬ 
niques.  Lastly,  with  regard  to  time-dependent  tension  failure  or  damage, 
additional  theoretical  and  experimental  work  is  needed.  A  tension  visco-damage 
model  recently  proposed  by  Whitman  (18)  shows  some  promise  in  this  area. 
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APPENDIX 


PROGRAM  VPDRVR:  INSTRUCTIONS /DOCUMENTATION 

This  Appendix  provides  input  instructions  and  documentation  for  the  VPDRVR 
computer  program  (FORTRAN  IV) .  VPDRVR  exercises  the  CAP75  viscoplastic  model 
developed  in  this  report  for  general  states  of  strain  or  stress  loading  schedules. 
The  general  solution  strategy  parallels  the  algorithm  presented  in  Table  2  for 
strain  loading  or  Table  3  for  stress  loading.  The  procedure  for  updating  the 
cap  hardening  parameters  follows  the  algorithm  presented  in  Table  4. 

Part  I  contains  user  input  instructions.  Part  II  describes  program  organi¬ 
zation,  subroutines,  and  program  variables.  Part  III  provides  input/output  for 
a  simple  benchmark  problem. 

Input  data  cards  are  grouped  in  the  following  categories: 

A.  (Cards  1  and  2) :  Heading  and  master  control 

3.  (Cards  3,4,  and  5):  Elastic  functions/parameter s 

C.  (Cards  6,7,3,  and  9):  Plastic  functions/parameters 

D.  (Card  10):  Viscous  functions/parameters 

Z.  (Cards  11,  12) :  Loading  schedules  for  stress  or  strain 
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Part  I.  USER  INPUT  INSTRUCTIONS 

A.  Problem  Initiation,  Heading  and  Master  Control  Cards. 
Card  1.  (15A4)  Heading 


Columns 

Variable 

Entry  Description 

Notes 

01-60 

(15A4) 

TITLE 

Descriptive  problem  title,  (program 
terminates  if  TITLE (1)  -  STOP) . 

(1) 

Card  2 . 

(415,  Al, 

2F10.0)  Master  Controls 

Columns 

Variable 

Entry  Description 

Notes 

01-05 

(15) 

LTYPE 

Loading  type  identification: 

”  0,  strain  loading. 

-  1,  stress  loading. 

(2) 

06-10 

(15) 

NTSEG 

Number  of  time  segments  to  define 
loading,  (Default  =  1,  Maximum  *  30). 

(3) 

11-15 

(15) 

ITMAX 

Number  of  Newton-Raphson  iterations, 
(Default  =»  10)  . 

(4) 

16-20 

(15) 

KPRINT 

Output  print  control:  (5) 

=  0,  standard  response  output. 

*  1,  above  plus  iteration  parameters 
■  2,  above  plus  yield  function  values. 

=  3,  above  plus  iterative  correction  vector 
^  4,  above  plus  Jacobian  matrix. 

20-21 

(Al) 

IPLOT 

Plot  control  for  response  data  written 
to  unit  11: 

=  Y,  (YES)  Data  written  to  unit  11 
*  N,  (NO)  Not  written 

(6) 

22-31 

(FI0.0) 

THETA 

Crank  Nicolson  integration  parameter,  9; 

o  <  e  <  i.o. 

(7) 

32-41 

(F10.0) 

CONVRG 

Convergence  tolerance  for  Newton-Raphson 
iteration,  (Default  =  0.01,  i.e.  1% 
relative  error) . 

(8) 
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B.  Elastic  Function  and  Parameter  Cards 


Card  3 . 

(215)  Selection 

of  Elastic  Functions 

Columns 

Variable 

Entry  Description 

01-05 

(15) 

IFBMOD 

Selection  of  bulk  modulus  function,  K(J.) : 

■  1,  K( J. )  =  BDATA(l) ,  (linear). 

-  2,  K(jp  -  BDATA(  1)  /  (1  -BDATA(2)  )  * 

(1-BDATA(2  )  *EXP  (BDATA(  3)  *J1) 

(Default  =*  1) 

06-10 

(15) 

IFSMOD 

Selection  of  shear  modulus  function,  G(J2) : 
=•  1,  G(J2)  =  SDATA(l) ,  (linear) 

*  2,  G(J2)  =  SDATA(l) / (1-SDATA(2)) * 

( 1-S DATA ( 2 ) *EXP ( -S DATA ( 3 ) * J2 ) ) 

(Default  =  1) 

Card  4. 

(7F10.0) 

Bulk  modulus  parameters,  BDATA. 

Columns 

Variable 

Entry  Description 

01-10 

(F10.0) 

BDATA(l) 

First  bulk  modulus  parameter. 

11-20 

(F10.0) 

3DATA(2) 

Second  bulk  modulus  parameter. 

21-30 
( F10 . 0) 

BDATA( 3) 

Third  bulk  modulus  parameter. 

Card  5 . 

(7F10.0) 

Shear 

modulus  parameters,  SDATA. 

Columns 

Variable 

Entry  Description 

01-10 

(F10.0) 

SDATA(l) 

First  shear  modulus  parameter. 

11-20 

(F10.0) 

SDATA(2 ) 

Second  shear  modulus  parameter. 

21-30 

(F10.0) 

SDATA(3) 

Third  shear  modulus  parameter. 

Notes 

(9) 


(10) 


Notes 

(11) 


Notes 

(12) 
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Plascic  Function  and  Parameter  Cards 


Card  6. 

Columns 

01-05 

(15) 


06-10 

(15) 


11-15 

(15) 


16-20 

(15) 


21-30 

(G10.0) 


Card  7. 

Columns 

01-10 

(F10.0) 

11-20 

(F10.0) 

21-30 
(FIG. 0) 


(415,  G10.0)  Selection  of  CAP75  functions 


Variable  Entry  Description 


IFFAIL 


Selection  of  failure  surface  function: 

fr" 


-  1,  g  =  -FDATA(l)  +  FDATA(2)*J1. 

F1 

-  2,  g  -  -FDATA(l)  +  FDATA(2) * 

¥1  EXP(FDATA(3) *J1) . 

(Default  =  1) 

IFCAPR  Selection  of  cap  surface  ellipse  ratio  R: 

■  0,  No  cap,  just  failure  surface. 

-  1,  R  -  CDATA(l) . 

-  2,  R  =*  CDATA(l)/  (1  +  CDATA(2) )  * 

(1.0  +  CDATA ( 2 ) *EXP ( CDATA ( 3 ) *EL) ) 

IFHARD  Control  of  cap  hardening: 

m  0,  No  hardening,  stationary  cap. 

=  1,  CAP75  hardening  function  is  used: 
e  =  W*(EXP(D*X)  -  1)  . 

W  *  HDATA(l) 

D  *  HDATA(2) 

KAPTYP  Selection  for  soil  or  rock  hardening  laws: 

*  0,  soil  material. 

*  1,  rock  material. 

XINITL  Initial  locatior  of  cap  X  on  axis. 


(7F10.0)  Failure  Surface  Parameters,  FDATA. 
Variable  Entry  Description 

FDATA(l)  First  failure  surface  parameter. 

FDATA(2)  Second  failure  surface  parameter. 

FDATA(3)  Third  failure  surface  parameter. 


Notes 

(13) 


(14) 

(15) 

(16) 

(17) 

Notes 

(18) 
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*  Card  8.  (7F10.0) 

Cap  Surface  Parameters  for  R,  CDATA. 

Columns 

Variable 

Entry  Description 

Notes 

01-10 

(F10.0) 

CDATA(l) 

First  cap  R  parameter. 

(19) 

11-20 

(F10.0) 

CDATA(2) 

Second  cap  R  parameter. 

21-30 

(F10.0) 

CDATA(3) 

Third  cap  R  parameter. 

*  Card  9.  (7F10.0) 

Hardening  cap  parameters,  HDATA. 

Columns 

Variable 

Entry  Description 

Notes 

01-10 

(F10.0) 

HDATA(l) 

First  hardening  paramter,  W. 

(20) 

11-20 

(F10.0) 

HDATA( 2 ) 

Second  hardening  parameter,  D. 

★ 

Skip  Cards  8 

and  9  if 

IFCAPR  *  0. 

D.  Viscous 

Function  and  Parameter  Card 

Card  10. 

(15,  3F10.0)  Selection  of  viscous  function/ parameters 

Co  Lunins 

Variable 

Entry  Description 

Notes 

01-05 

(15) 

IFVISC 

Selection  of  viscous  function  r> : 

-  1,  $  -  ( f  / ANORM)  **EXPN . 

-  2,  $  »  EXP (([/ ANORM) **EXPN)-  1. 

(Default  =*  1) 

(21) 

05-15 
(FI 0.0) 

EXPN 

Exponent  in  ?  function, 

(Default  -  1.0)  . 

(22) 

16-25 

(F10.0) 

GAMMA 

Fluidity  parameter,  y. 

(23) 

26-35 

(F10.0) 

ANORM 

Normalizing  constant  in  j  function, 
(Default  =  max(FDATA(l) ,  0.01) 

(24) 
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E.  Input  Loading  Schedule  and  Time  Steps. 

Repeat  card  set  11  and  12  NTSEG  times;  NS  ■  1,  NTSEG 


Card  11 

(F10.0,  215) 

Time  segment,  number  of  stops,  print  control. 

Columns 

Variable 

Entry  Description 

Notes 

01-10 

(F10.0) 

TS(NS) 

Time  at  end  of  segment  NS. 

(25) 

11-15 

(15) 

NUMDT(NS) 

Number  of  times  steps  within  time 
segment  NS . 

(Default  =  10) 

(26) 

16-20 

(15) 

IPRNT(NS) 

Print  interval  for  standard  output: 

=  1,  every  time  step  prints  output. 

*  n,  every  nth  step  prints. 

(Default  «  1) 

(27) 

Card  12 

(6F10.0)  Stress  or  strain  load  vector  at  time  TS(NS). 

Columns 

Variable 

Entry  Description 

Notes 

01-10 

(F10.0) 

PLOAD(l.NS) 

*11  (or  en)  at  TS(NS). 

(28) 

11-20 

(F10.0) 

P LOAD (2 , NS) 

322  ^°r  Zn)  3t  TS^NS^ • 

21-30 

(F10.0) 

? LOAD (3, NS) 

a33  (or  e33)  at  TS(NS) . 

31-40 
(F10. 0) 

P LOAD (4, NS) 

*1,  (or  e  )  at  TS (NS) . 

41-50 

(F10.0) 

PLOAD( 5, NS) 

J13  (or  e13)  aC  TS(NS)- 

51-60 
( F10 . 0) 

P LOAD (6, NS) 

J,,  (or  c  )  at  TS (NS)  . 

Z-i  23 

***END 

OF  INPUT  FOR  ONE  PROBLEM*** 
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Commentary  Notes  with  Input  Instructions: 

1.  Problems  may  be  run  back-to-back.  Terminate  the  last  problem  by  writing 
STOP  in  columns  1  to  4. 

2.  Strain  loading  implies  the  six  components  of  strain  will  be  specified 
individually  during  the  loading  schedule.  Similarly,  stress  loading 
implies  the  six  components  of  stress  will  be  individually  specified. 

3.  For  either  stress  or  strain  loading,  NTSEG  is  the  desired  number  of 
time  segments  to  define  the  loading  histories  in  a  piecewise  linear 
fashion. 

4.  Generally  10  iterations  is  more  than  sufficient  to  achieve  convergence. 

If  convergence  is  not  achieved,  it  is  a  strong  indication  that  the  time 
step  is  too  large.  Note  that  convergence  of  the  Newton-Raphson  procedure 
does  not  guarantee  accuracy.  Accuracy  can  only  be  assured  by  repeatable 
solutions  with  smaller  time  steps. 

5.  Standard  output  includes  stress  or  strain  responses,  cap  location,  number 
iterations  to  converge,  stress  invariants,  and  type  of  response.  For 
KPRINT  >  0,  additional  information  is  given  primarily  for  debugging  purposes. 

6.  Standard  response  data  is  written  to  unit  11  for  subsequent  plotting  on  a 
CALCOMP  plotter.  Subroutine  GRAPH  is  used  for  plotting  and  maybe  removed 
or  replaced  if  desired. 

7.  For  THETA  =  0.0,  the  solution  algorithm  is  explicit  resulting  in  linear 
equations  (i.e.  no  Newton-Raphson  iteration).  For  THETA  >  0,  the  algorithm 
is  implicit  and  generally  more  accurate  for  a  given  time  size,  but  requires 
Newton-Raphson  iteration.  For  THETA  _>  0.5,  the  algorithm  is  unconditionally 
stable. 

8.  The  convergence  tolerance,  CONVRG,  is  tested  against  the  ratio  formed  by 
the  norm  of  the  correction  vector  for  stress  (or  strain)  divided  by  the 
norm  of  the  stress  (or  strain)  vector.  Norms  are  Euclidean. 

9.  The  nonlinear  bulk  modulus  function  given  by  IFBMOD  =  2  is  taken  from 
CAPDRIVER  (Reference  17).  It  is  a  function  of  Ji  (first  stress  invariant) 
and  is  treated  the  same  for  loading  or  unloading.  Additional  functions 
may  be  added  to  program  in  FUNCTION  DI(I,J). 

10.  The  nonlinear  shear  modulus  function  given  by  IFSMOD  =  2  is  a  function  of 
J2»  second  deviator  stress  invariant  (see  Note  9) . 

11.  For  future  program  expansion,  3DATA  is  dimensioned  to  7  to  allow  incorpora¬ 
tion  of  higher  order  nonlinear  functions. 

12.  SDATA  is  dimensioned  to  7  (see  above). 

13.  For  IFF AIL  *  1,  the  failure  surface  is  standard  Drucker-Prager  (or  Von 
Mices  if  FDATA(2)  «  0.0).  For  IFFAIL  =  2,  the  failure  surface  is  the 
exponential  form  suggested  by  Sandler  for  CAP75.  Additional  functional 
forms  may  be  added  to  the  program  in  FUNCTION  FG1. 
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14.  By  setting  IFCAPR  -  0,  the  plasticity  model  is  governed  by  only  the  failure 
surface.  For  IFCAPR  *  1  or  2  the  cap  surface  is  included  with  R  given  by 
the  corresponding  functional  form.  Additional  functional  forms  for  R  may 
be  added  to  program  in  FUNCTION  FRCAP.  (Note  for  IFCAPR  «  2,  R  m  R(EL) 
where  EL  is  "L"  of  cap) . 

15.  If  desired,  a  nonhardening  cap  surface  may  be  used  by  setting  IFHARD  *  0. 
Otherwise  the  CAP75  hardening  function  is  employed-  New  hardening  functions 
can  be  employed  by  modifying  SUBROUTINE  CAP75. 

16.  See  Table  4  for  the  special  hardening  rules  for  soils  (KAPTYP  *  0) . 

17.  The  initial  location  of  X  defines  the  starting  position  of  the  cap  surface. 

The  program  checks  that  XINITL  is  not  greater  than  FCUT,  i.e.  the  inter¬ 
section  of  the  failure  surface  with  Jl  axis.  If  it  is,  XINITL  is  auto¬ 
matically  reset  slightly  less  than  FCUT.  Note,  the  so-called  Von  Mises 
Transition  employed  by  Sandler  is  not  included  in  this  development.  Thus, 
if  it  is  desired  to  obtain  steady-state  viscoplastic  solutions  to  exactly 
match  CAP75  plasticity  solutions,  XINITL  should  be  chosen  so  that  the 
initial  L  location  is  not  greater  than  zero. 

18.  The  "standard  Sandler"  CAP75  failure  surface  is  the  form  given  by  IFFAIL  * 

2.  In  which  case  FDATA(l)  =*  A,  FDATA(2)  =  C,  and  FDATA(3)  =*  B. 

19.  The  "standard  Sandler"  CAP75  cap  surface  parameter  is  the  form  given  by 
IFCAPR  =  1,  i.e.,  CDATA(l)  =  R. 

20.  If  IFHARD  =  0,  HDATA(l)  and  HDATA(2)  are  read  but  not  used.  If  IFCAPR  - 

0,  cards  S  and  9  are  not  read.  HDATA  as  well  as  FDATA  and  CDATA  are 
dimensioned  to  7  for  future  program  expansion.  » 

21.  For  geological  materials  IFVISC  =  1  is  generally  the  most  popular  form 
for  the  viscous  function.  Additional  functional  forms  such  as  Equations 
60  and  61  may  be  added  to  the  program  in  SUBROUTINE  PHIF. 

22.  EXPN  need  not  be  a  whole  number,  but  must  be  greater  than  zero. 

23.  GAMMA  has  units  of  inverse  time,  the  units  (e.g.  seconds,  hours,  years) 

correspond  to  the  loading  time  units  TS  in  Card  11. 

24.  Generally  the  default  value  of  ANORM  is  appropriate  providing  FDATA(l)  +  0.0. 

AN0RM  should  not  be  viewed  as  an  independent  material  parameter  since  it  is 
always  associated  with  GAMMA  in  the  quotient  GAMMA/ ANORM**EXPN . 

25.  Up  to  30  time  segments  may  be  used  to  define  a  piecewise  continuous  collection 
of  straight  lines  to  define  loading.  For  the  first  time  segment,  the  program 
automatically  assumes  initial  time  is  zero,  i.e.  TS(0)  =  0.0.  Thus,  TS(1) 

is  the  time  at  the  end  of  first  segment,  TS(2)  is  the  time  at  the  end  of 
the  second  segment,  etc.  Successive  values  of  TS(NS)  must  be  greater  than 
the  previous  value. 
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26.  Any  number  of  time  steps  may  be  assigned  to  each  time  segment.  Accuracy/ 
stability  is  controlled  by  the  time  step  size  so  that  it  is  good  practice 
to  repeat  solutions  by  doubling  the  value  of  NTJMDT(NS)  .  Although  the 
time  step  size  may  be  specified  differently  in  each  time  segment,  it  is 
good  practice  not  to  make  changes  in  At  between  segments  by  a  factor  of 
more  than  2. 

27.  The  printout  interval  may  be  specified  differently  for  each  time  segment. 

28.  Loading  values  at  the  end  of  each  time  segment  are  specified  individually 
for  each  vector  component  of  strain  if  LTYPE  *  0,  or  each  vector  component 
of  stress  if  LTYPE  »  1.  For  the  first  time  segment  the  initial  loading 

and  responses  are  automatically  assumed  zero  i.e.,  a(0)  *  e(0)  ■  0.  Standard 
continuum  mechanics  sign  convenctions  are  observed  for  all  input  and  output. 
For  example,  if  a  uniaxial  stress  loading  cycle  is  desired  in  which  cr^  is 
compressed  at  a  constant  rate  to  a  stress  value  -10.0,  held  constant,  then 
reverse  loaded  at  a  constant  rate  to  a  tensile  stress  value  of  +1.0,  and 
again  held  constant;  we  infer  NTSEG  =  4,  and 


a  is  described  by: 


PL0AD(1,1) 
PL0AD(1, 2) 
PL0AD(1,3) 
PL0AD(1, 4) 


=■  -10.0 
-  -10.0 
=  +1.0 
*  +1.0 


and  all  other  stress  components  (PL0AD)  are  zero. 


PART  II.  VPDRVR  Documentation 


Table  Al  describes  subroutines  and  functions  employed  in  VPDRVR  along  with 
associated  calls.  Table  A2  illustrates  the  program  flow  path. 

Listed  below  are  all  COMMON  statements  with  a  description  of  their  variables 


1.  COMMON/MASTER/ 

THETA  ■  9,  Crank-Nicolson  integration  parameter 

CONVRG  =*  convergence  tolerance  ratio 

DT  =  At,  current  time  step  size 

LTYPE  =  loading  type;  strain  -  0,  stress  =  1. 

NTSEG  *  number  of  time  segments 
KPRINT  =  print  control  parameter 

ITMAX  *  maximum  number  of  Newton-Raphson  iterations 

2 .  COMMON /ELAST/ 

BDATA(7)  =  bulk  modulus  coefficients 
SDATA(7)  =  shear  modulus  coefficients 
IF3MOD  =  functional  form  of  bulk  modulus 
IPS MOD  =  functional  form  of  shear  modulus 

3.  COMMON/PLAST/ 

X  *  cap  surface  location  X 
EL  *  cap  surface  location  L 

VH  »  £,  viscoplastic  volumetric  strain  hardening 
VHMAX  =  Ejnax’  saaximum  value  of  e 
FDATAi.7)  =  CAP75  failure  surface  constants 
CDATA(7)  =  CAP75  cap  surface  constants 
HDATA(7)  =  CAP75  cap  hardening  constants 
IFFAIL  ■  functional  form  of  failure  surface 
IFCAPR  =  functional  form  of  cap  surface 
IFHARD  ■  functional  form  of  cap  hardening 

KAPTYP  ■  indicator  for  rock  or  soil  hardening  law;  soil  =  0,  rock  *  1 

4.  COMMON/ VI3C0/  1 

EXPN  =»  N,  exponent  in  j  function 

GAMMA  *  v,  fluidity  parameter 

ANORM  *  fg,  normalizing  constant  in  p  function 

IFVTSC  =  functional  form  of  p  function 

5.  COMMON/ RESULT/ 

EPS (6)  =  z,  strain  vector 
SIG(6)  »  cr,  stress  vector 
EVP (6)  *  s  ,  viscoplastic  strain  vector 

-vp 
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EVPD0T(6)  -  e  ,  viscoplastic  strain  rate  vector 
~vp 

SJ1  -  J^,  first  stress  invariant 

SJ2  ■  J 2»  second  deviator  stress  variant 

ITER  *  iteration  number 

6.  COMMON / LOADN G 

TS(30)  =  time  at  end  of  time  segment 
PL0AD(30)  =  loading  at  end  of  time  segment 
IPRINT(30)  *  print  interval  in  time  segment 
NUMDT(30)  ■  number  of  time  steps  in  time  segment 

7.  COMMON / UTLTY / 

P(6)  *  P  or  P  in  Equation  72  or  79 
PP(6,6)  «  ,  Jacobian  matrix  (Equation  101) 

RHS(6)  *  q  or  q  in  Equation  73  or  78 

CELC0R(6)*  6o  or  5e,  correction  vector  for  Newton-Raphson 
DPL0AD(6)  ■  A q  or  Ae,  load  increment 

DEVPKK  »  Aoi,  increment  of  volumetric  viscoplastic  strain 

8.  COMMON/YIELDF/ 

G1  -  gj_,  value  of  plastic  surface  function 
G1P  -  g^,  first  derivative  of  w.r.t.  J.^ 

G1PP  -  g!^,  second  derivative  of  w.r.t. 

G2  ■  g2>  value  of  plastic  surface  function  g2 
G2P  ■  g2>  first  derivative  of  g2  w.r.t.  J., 

G2PP  ■  g'2>  second  derivative  of  g2  w.r.t.  J2 

9 .  COMMON /VFLOWF/ 

F  ■  f  or  fc>  yield  function  value 

PHI  *  viscous  function  value 

PHIP  -  $' ,  derivative  of  <J>  w.r.t.  f 

PHIPRE  ■  $n,  value  of  $  at  end  of  previous  time  step 

VHPRE  «  tn,  value  of  £  at  end  of  previous  time  step 

ISURF  *  indicator  for  governing  surface:  failure  surface  =  0,  cap  =  1 
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Table  Al.  Subroutines  and  Functions 


Main  program:  Executive  duties  for 
controlling  input  calls,  initialization, 
time  step  sequencing,  and  output. 


Reads  in  functional  forms  and  para-  VPDRVR 

meters  for  elastic,  plastic,  and  viscous 
components  of  viscoplastic  model. 

Reads  in  loading  schedule,  time  step,  VPDRVR 
and  print  controls. 

Determines  stresses  for  strain  loading  VPDRVR 

by  Newton-Raphson  iteration  and  updates 
responses. 

I  Determines  strains  for  stress  loading  VPDRVR 

by  Newton-Raphson  iteration  and  updates 
responses . 

Computes  stress  invariants,  yield  GETSIG 

function  value,  viscoplastic  strain*  GETEPS 

i  rate,  and  forms  Jacobian  for  GETSIG. 

i  ; 

!  Computes  yield  functions  gi  and  g2  VPLAST 

'  and  their  derivatives  for  CAP75 
|  plasticity  and  updates  hardening 

j  parameters. 

Computes  viscous  flow  function  $  and  VPLAST 

,  its  derivative  .  ! 

!  ! 

j  Gauss  elimination  equation  solver  GETSIG 

Calcomp  plotting  subroutine  >  VPDRVR 

I 

Determines  components  of  D  elastic  !  GETSIG 

matrix  I  GETEPS 


Computes  g. ,  g'  and  g"  as  a  function  j  NPUT4M 

, .  1  .  1  I  CAP75 

first  stress  invariant.  j 

Computes  cap  parameter  R  as  a  function  j  CAP75 

of  cap  location  !  NPUT4M 


NPUT4M 

LOADS 

GETSIG 

GETEPS 

GRAPH 

FG1 

FRCAP 


SOLVER 

VPLAST 

DI 

VPLAST 

DI 


CAP  7  5 
PHIF 


FG1 

FRCAP 


Table  A 2.  VPDRVR  Program  Flow 


[  PROGRAM  VPDRVR 

1 

1  SUBROUTINES 

j  (main) 

1  FUNCTIONS 

1 

(a)  Master  Control 

(b)  Material  Properties  * - 

(c)  Loading  Schedule  •* - 

INITIALIZE  VARIABLES 

TIME  LOOP 

(a)  Determine  load  increment 

(b)  Get  solution  +  update  " 

(c)  Print  results 


Return  to  INPUT  for  new 
problem. 


NPUT4M 

LOADS 


GETS I G 
(or) 
GETEPS 


FG1  I 
FRCAP! 


SOLVER 


V? LAST  1 


CAP75 

-J 

n 

FRCAP 

PHIF 
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PART  III  Example  input/output 


The  uniaxial  strain  problem  presented  in  the  main  body  of  this  report  is 
used  here  as  a  benchmark  example  for  the  VPDRVR  program.  On  the  following 
page  is  listed  the  formatted  input  (card  images)  corresponding  to  the  input 
instructions  in  Part  I  of  this  appendix.  Subsequent  pages  show  the  output 
from  the  VPDRVR  program  including  a  restatement  of  the  input  with  default 
values,  as  well  as,  response  data.  For  brevity,  the  printout  of  response 
data  is  limited  for  each  of  the  four  time  segments  to  include  only  the  time 
segment  midpoint  and  end  values.  All  data  is  adequately  labeled  and  should 
cause  no  confusion  except,  perhaps,  for  the  nomenclature  associated  with 
the  response  "STATE".  Here  the  following  meanings  are  implied: 

CAP-VP  -  viscoplastic  flow  above  cap  surface 
CAP-SS  »  steady  state  response  on  (or  near)  cap  surface 
FSURF-VP  »  viscoplastic  flow  above  failure  surface 
FSURF-SS  -  steady  state  response  on  (or  near)  failure  surface 
ELASTIC  *  stress  state  is  in  elastic  domain 
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Input  to  VPDRVR,  Uniaxial  Strain,  y  =  0.01. 


00010 

SAMPLE  OUTPUT  OF 

VPDRVR  -  UNIAXIAL  STRAIN  LOAD 

00020 

0  4 

ON  .75 

00030 

1  1 

00040 

66.67 

00050 

40.0 

00060 

2  1 

1 

0  -.1888 

00070 

0.25 

0.18 

0.67 

00080 

2.5 

00090 

0.066 

0.67 

00100 

1 

0 

.01 

00110 

1  . 

80 

40 

00120 

-0.03 
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5  • 

320 

160 

00140 

-0.03 

00150 

tr  «r 

U  ♦  «J 

40 

20 

00160 

-0.0225 

00170 

7.5 

160 

80 

00180 

-0.0225 

00190 

STOP 
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Output  from  Vl’OKVK  (5  pages) 
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